Proteomic analysis reveals key differences between squamous cell carcinomas and adenocarcinomas across multiple tissues

Squamous cell carcinoma (SCC) and adenocarcinoma (AC) are two main histological subtypes of solid cancer; however, SCCs are derived from different organs with similar morphologies, and it is challenging to distinguish the origin of metastatic SCCs. Here we report a deep proteomic analysis of 333 SCCs of 17 organs and 69 ACs of 7 organs. Proteomic comparison between SCCs and ACs identifies distinguishable pivotal pathways and molecules in those pathways play consistent adverse or opposite prognostic roles in ACs and SCCs. A comparison between common and rare SCCs highlights lipid metabolism may reinforce the malignancy of rare SCCs. Proteomic clusters reveal anatomical features, and kinase-transcription factor networks indicate differential SCC characteristics, while immune subtyping reveals diverse tumor microenvironments across and within diagnoses and identified potential druggable targets. Furthermore, tumor-specific proteins provide candidates with differentially diagnostic values. This proteomics architecture represents a public resource for researchers seeking a better understanding of SCCs and ACs.

In the case of the overall survival/disease free survival there are cohorts collected at different times having different time to followup. Therefore, comparing outcomes may be difficult, as censoring may be primarily a late or early event, and may be influenced by date of collection (due to therapy). There does not appear to be multiple testing correction methods applied to survival p-values. Given the large number of exploratory tests for outcomes, multiple testing correction may be appropriate to consider.
Specific Items: Lines 125-126: It would be helpful in showing KRT5 and TP63 to again find a quantification of the assertion that they are "high and ubiquitous". There is some variation in KRT5, at least. It appears that many thyroid cancers have lower KRT5, so it might be helpful to summarize observations with quantification (perhaps by tissue of origin). For instance, how does the variability in AKT1 or KRT5 compare (for instance, what quantile of variability using MAD as done elsewhere)?
On the basis of the methods and the color gradient in the figure 1d, blue is "missing" data. If this is the case, it may be worth pointing this fact out. FGFR3, PTEN and TP53 may be mutated and thus either not sufficiently translated or translated with alternate forms. Given the (relatively) low expression of those proteins (FGFR3/PTEN) in samples, it may be worth exploring this further.
Identification of proteins may relate back to tissue specificity so it would be helpful to know if the total numbers were significantly different across tissues. Specifically on line 115 "no major differences in the coverage between the 17 SCCs". I would encourage a test of this assertion. Figure 1c suggests there are significant differences. Were samples run by tissue of origin? If they were randomized, then the differences in total protein could be biologically meaningful and worth further investigation. This is particularly true with copy number changes that have been seen in squamous cell cancers.
Batch correction: When batch correction was done, there was missing data present. Was the missing data imputed in the batch correction or left missing? That is, the missing data is imputed as a low expression however it is possible that batch correction could increase the level of expression of these missing values? This would presumably introduce an unintended batch effect in that case?
Commonly expressed proteins in SCC or AC: It is not clear to me, when 5130 commonly expressed SCC proteins and 4845 commonly expressed AC proteins yield 5838 common proteins of SCC or AC. Is this perhaps the union of the two sets? It would seem more common to assume that the intersection of the two is used as the basis for common proteins. Overall, there are two potential effects: 1) presence of proteins in one of AC/SCC and 2) expression differences among proteins in both AC/SCC. It is not clear which (or both) of these is being addressed. It is expected that AC and SCC have unique proteins (of interest in and of itself) but it is not clear how batch correction for missing proteins can be successfully done unless it is implicitly a form of imputation.
Survival Analysis: Gene expression was dichotomized using the maxstat approach. There is no multiple testing corrections performed when testing targets in multiple diseases. For example, Fig 2g. Given the large numbers of tests and calling out only specific targets in specific diseases, it is difficult to assess the statistical significance of these findings. Further, there are references for the use of maxstat that should be considered, vs a URL: e.g.,  : This applies to the PLIN1 finding, but likely applies to a number of other findings. In the case of PLIN1, it appears that the large log2 ratio for rare SCC comes from the fact that most (188/228) common samples did not identify PLIN1. In these cases, it would be helpful to provide any information on presence of peptides associated with PLIN1 in the common cancers. It would be extremely interesting if this protein was only (or substantially) present in rare squamous tumors, potentially due to a deletion in "common" squamous tumors. However due to label-free characteristics it is possible that technical artifacts caused this difference (perhaps fewer proteins overall, run order, peptide interference).  Figure 3: This is difficult to understand. For instance, it appears that SPRR1A is high in normal esophagus from the human protein analysis, yet missing (it appears) in the pan-squamous cohort. Please clarify this finding.
Supplemental Figure 5c: There are several tissue types labeled with color (blue, red). Does this color indicate a particular relationship? It does in the case of Thyroid and Gallbladder (similar) but it's not clear with breast scc. Line 160,161: It is not clear if the term overrepresented is accurate. It appears to be the results of a wilcoxon test, therefore it is overexpressed or overabundant. : It is not clear what data is used to assert that prognosis "were significantly affected by these DEPs". : This appears to be a statement about other work, not specifically results from this cohort although it is not clear.

-Characterization of HPV-related SCCs
Lines 402-403: it seems likely that the positive rate above 54% is in the HPV+ patient cohort (as is HPV18 6.67), but it is not clear.
Lines 410-411: p53 loss occurred in HPV-cases as well. This would be an example where a statistic may provide more information (perhaps the frequency of p53 loss).
Line 411: "CDKN2A expression was highly correlated" -please include p value or correlation coefficient. Figure 6b does not include these numbers either.
Line 418: "We identified 8 patterns of differential pathway ..." -was tissue of origin included in the limma analysis for this? Was this driven by tissue of origin?  Clinical Sample Acquisition (582). Patients were excluded if they had advanced disease (perhaps this could include clinical stage for clarity). Also, "any condition that may influence the outcome evaluation". It would be helpful to indicate that some patients are missing outcomes. Presumably these were dropped from survival analysis, but may have been censored at time 0. It would be helpful to clarify how this was addressed.
Additionally, samples were collected from 2001 to present. The manuscript states "all patients provided written consent". This suggests that a general banking protocol was used to collect tissues, distinct from this specific study. It would be helpful to clarify this. In the case of the adenocarcinomas, it is explicity stated that patients were consented prior to surgery.
Patients with adenocarcinomas were excluded if radiation was used (unclear if this was preoperative). Presumably this was not the case in the squamous cohort, was this evaluated?
Database Searching (640): What version of Uniprot was used? It may be present in a figure, but not stated in the methods. Also, were missed cleavages considered for trypsin identification? Modifications? More details here would be needed to understand how the searching was done.
MAD-based protein selection: line 759-767: What does "the 20% bottom mad proteins of each organ were combined and duplicates removed" mean? It would seem reasonable to remove the bottom 20% of mad proteins due to lack of variability, but it's not clear what the term "combined" means here (averaged?).
line 785-794.: Why would the 1% bottom mad proteins be selected? These are the least variable (assuming smallest mad values are meant). It is not clear why the 1220 proteins were used to compute an average expression value (for what purpose). In particular, ordering by mad does not impact the random forest. Perhaps something different was meant, but it was not clear from the description. In general, it would also be good to indicate if all filtering and protein selection was done only on the training set (independent evaluation) or on both the training and testing sets. line 721: "multiple testing correction using ..." Line 811/812: Please define HPV16 as the "main type".
Discussion comments: line 511: "involved in rare SCC initiation". It would be more accurate to describe them as regulators involved in rare SCC since it is a descriptive finding.

Reviewer #2 (Remarks to the Author): Expert in SCC proteomics
In the manuscript entitled "Proteomic Landscape of Pan-squamous Cell Carcinomas", Song et al. provided an interesting panorama of SCC proteomes from 333 patients and 17 sites using mass spectrometry label-free quantitation. By performing several group comparisons (SCC vs AC; common SCC vs rare SCC; HPV+ anogenital SCC vs HPV-anogenital SCC), the authors indicated signatures of differentially abundant proteins that have a prognostic significance and are able to modulate specific pathways. Four groups of SCC types were defined by clustering tumor proteomes and six subtypes of SCCs were determined based on immune enrichment of cell composition. Molecules from these groups may potentially modify biological pathways and are promising as druggable targets. Finally, a machine learning model was built to predict the tumor site of origin based on SCC proteomes and 3 proteins were selected for validation using IHC.
Overall, that is a detailed study that provided insights into the biology of a common cancer type (SCCs) by profiling the tumor proteomes followed by an extensive bioinformatic analysis. The significant number of patients evaluated is certainly a highlight of the study. Since large-scale analysis of tumors available in the literature have focused on DNA/RNA levels, describing the proteome and the potential functional implication is definitely novel and opens up new perspectives for the study of cancer. The description of a protein signature that can possibly define the origin of SCCs is also remarkable and reveled targets with potential to translate to the clinics. For instance, the authors provided a valuable resource to the scientific communities for further exploration. Considering the methodology, a noteworthy point is that the manuscript presents a deep characterization of SCCs achieved by implementing a strong and robust bioinformatic analysis. The proteomic results indicate a successful workflow of sample preparation, MS run and data analysis (~ 15,000 protein groups quantified for SCCs and ~10,000 for ACs), even though some additional information should be provided to assure the reproducibility and quality of the data.
In summary, the study is of major interest for the oncology and proteomics fields and is in-line with most articles published in Nature Communications. There are therefore some shortcomings that should be addressed.

Results
1. Line 109: The paper relies on proteomics data from a large group of SCC samples and the quality control of MS runs is certainly one of the main concerns to assure the reliability of bioinformatic analysis. However, I am not sure about QC analysis. How was the set of repeated samples selected for correlation (SFig. 2b)? Additionally, the criteria for QC of individual samples were not described in the text and it would be good to see some QC results for all MS runs. Plotting Spearman's correlation, counting valid values per sample or, if suitable, describing retention times or m/z for trypsin autolysis peaks across samples are suggested approaches.
2. Normalized protein intensities from pan-SCC and pan-AC were log2 transformed. Even that is not always true, log transformation reduces the skewness of large-scale data and make it more closely to a normal distribution. If that is the case, it would be appropriate to use parametric tests for group comparison instead of the non-parametric analysis described in the manuscript (Wilcoxon, Kruskal-Wallis). Did the authors test the normality assumption of the data? Please comment. 5. Lines 179 and 325: Besides using their own proteomics data for Cox analysis of DEPs between SCCs and ACs, 9 RNA databases from cancer tissues were used to assume the role of proteins in prognosis. Did the authors evaluate if transcripts are somehow correlated with proteomics data? Otherwise, the assumptions may not be true. 6. Fig. 5a: The authors should amend Subtype 2 name to FaSq. 7. Fig. 5d: Not mentioned in the main text. 8. Fig. 6b: How can the authors explain the high levels of protein pRB in HPV-infected tumors? Since the silencing of pRb by E7 viral protein produces a rise in p16, the abundance of pRB is not in agreement with what is reported in the literature or in this study. Also, SOX2 is not frequently associated with HPV infection and it would be appropriate to describe the rationale for including this protein in the analysis. 9. Fig. 6e: The error bars are too large and it is difficult to believe that there is a real difference between HPV+ and HPV-cases for these proteins, even with an adjusted p value. Maybe including the protein abundances for the non-anogenital HPV-negative SCC types would make the hypothesis of Fig.6f stronger.
10. Line 459: I didn't get the point of why the authors evaluated the 3 markers by IHC. If a signature of 19 proteins was accurately able to discriminate SCC tumor based on their origin (buy the way, this information is not stated in the main text), why the 3 proteins were selected for validation? I understand the validation phase is an important step in large-scale experiments, but analyzing the 3 proteins alone does not make sense in the context of the classifier and did not make the proteomic data stronger.
Discussion 11. It is not necessary to extensively re-state the key findings in the discussion.

12.
Line 580: References should be provided for WHO classification and TNM system. 13. Fundamental information is missing in the proteomic workflow and it is difficult to judge the reproducibility and quality of the methods employed. Please provide additional information.
Sample preparation a. Why did the authors add an acetone precipitation step in the FASP protocol? FASP did not take care of contaminant removal? The authors should provide the appropriate references.
b. How was alkylation and quenching performed? c. The authors declare that a trypsin-to-protein ratio of 1:50 was used for digestion, but there is no information on protein quantification. How were protein levels determined? The authors also state that "Target on-column load was 200ng total peptide per injection", indicating that peptides were also quantified. Please clarify. d. It is appropriate to present specifications of the trypsin used for digestion. e. What is MS water? How were peptides acidified to stop digestion? LC-MS/MS f. Some fundamental aspects of the MS runs are missing, like m/z range, mode of data acquisition (DDA? DIA?), resolution, etc.
g. Although the identifier is provided in the text, the repository where raw files are deposited is not informed and data could not be accessed.

MS search
h. Details for the Uniprot library should be provided, including download date, number of residues considered. Were SwissProt and TrEMBL entries considered? i. A fragment ion tolerance of 0.05Da was used. Why did the authors use such a restrictive cut-off? I am afraid that important protein identification was lost.
j. How did the authors handle contaminants? k. I could not find any information about variable and fixed modifications, or the retention time window considered.
Data analysis l. Replacing missing values engenders intense debate in the scientific community. Are there any reasons why missing values were replaced? Maybe keep the original data would be appropriate and statistics would take care. Have other approaches been tested for data imputation to assure that replacing by one-tenth of the minimum intensity is the most suitable strategy?
14. I was wondering whether the separation of SCCs and ACs in PCA before batch effect correction just reflects their distinct biological characteristics. How can the authors be sure that the separation was a batch effect?
15. Not clear if IHC was performed in the same cohort as proteomics.
16. For HPV grouping, how the authors defined if HPV16 is the main type (group 2) or not (group 3) in multiple infections?
Understanding the molecular pathways driving histologically similar cancers across anatomic sites may provide new treatment paradigms that historically have been site-specific. Differences in mutational patterns and gene expression profiles between squamous cell carcinomas (SCC) and adenocarcinomas (AC) arising across anatomic sites have been well described using common resources such as TCGA. However, the proteomics landscape of SCC across anatomic sites has not been previously investigated in large numbers of tumors. The current manuscript describes proteomic patterns in 333 treatment-naïve squamous cell carcinoma (SCC) tissues obtained from 17 organ sites and 69 treatment-naïve adenocarcinoma (AC) tissues obtained from 7 organ sites from a single university-based hospital in Shanghai, China. Proteomic characterization of tissues was conducted using a mass spectrometry-based approach validated using tissue microarray (TMA) immunohistochemistry (IHC).
Major findings include the elucidation of pathways differentiating SCC from AC (keratinization, glucose metabolism and extracellular matrix), molecules within those pathways associated with disease prognosis, and proteomic clusters/immune subtypes that may represent potential druggable targets. The resulting data repository will serve as a unique and valuable shared resource for investigators to use in the future.
Methods and results are described in great detail, yet there are some key points that should be clarified and/or expanded upon.
Case selection and classification: • The methods state that the cases were randomly selected. How was this accomplished, and what percentage of the total SCC cases treated in the 18-year range do the cases included in the current study represent? Exclusion criteria are presented, yet it's unclear how many patients were excluded for the reasons listed.
• While the overall number of tumors (n=333) is substantial, the numbers of samples available per anatomic site ranged from 10-22 for SCC and 8-12 for AC. Therefore, inferences drawn for specific sites are limited by small sample size. This point should be added to the discussion.
• The classification of rare versus common tumors is unclear. The authors state that the WHO Classification of Tumors was used, yet there is no reference, and some cancers seem to be misclassified. For example, SCC of the vagina is very rare (i.e. incidence is less than 6 per 100,000), yet it is included here as a common cancer.
• As the authors point out in the Background, metastatic SCC's (or primaries with elevated metastatic potential) are an important clinical challenge. However, only primary SCC's were included in this case series, and no information was provided on whether or not patients developed metastases during follow-up. This is an important design limitation that should be discussed.
Patient follow-up and survival analysis: • No information is provided on the average length of follow-up (and range), as well as whether or not patients were lost to follow-up, and if so, how they were handled in the analysis. It is also not clear that the proportional hazards assumption was assessed.
• Why were patient age and gender (as well as other patient characteristics such as stage at diagnosis) not considered as potential covariates in the multivariable modeling, along with the three covariates stated (protein expression, organ and histology)? Results: • In general, it was difficult to follow the results section given the sheer number of figures, figure panels and supplementary materials. The figures were very detailed, as were the supplementary data (often patient-level data files). In many cases, it would have been helpful to create summary tables that allow the reader to directly compare percentages between groups and better ascertain the statistical significance of the observed results. For example, for Figure 1b-it would be helpful to show the information in tabular form so that percentages of samples across anatomic sites could be more directly compared with respect to tissue characteristics such as stromal score and keratinization; statistical significant testing could be used to determine which differences are most likely to be real and not due to chance. The raw data are included in supplementary Table 1, but a table showing the percentages across groups would be most helpful to view.
• Regarding the HPV results, it is difficult to glean from Figure 6e whether the protein patterns depicted are specific to HPV 16. It would be helpful to present HPV type-specific prevalence by tumor type in tabular form, and then present the percentages of each of the five HPV groups defined in Fig 6c that express proteins corresponding to the different pathway groups of interest defined in 6d. Were HPV16 E6 and E7 proteins detected in any of the SCC samples?
• The survival analysis described in Fig 2e is intriguing, given that these proteomic features may be useful prognostic markers. It appears as if fewer proteins were predictive of survival in the current PanSCC dataset compared to a majority of the TCGA datasets included in Fig 2e. Could this be a function of sample size? It would be helpful if the Hazard Ratios and 95% confidence intervals were provided to better interpret the magnitude and precision of these estimates.

Reviewer #1 (Remarks to the Author): Expert in bioinformatics
This paper is an important work that highlights the generation of a very large dataset consisting of 333 patients with squamous cell tumors from 17 different organs. The authors use label-free, mass spectrometry-based proteomics to characterize this cohort and compare with a similarly generated cohort of 69 patients with adenocarcinomas from 7 different organs. The generation of this dataset is a noteworthy and significant accomplishment as a mechanism for comparing the proteomics of squamous tumors from multiple organs consistently (same methodology, same instrument). In addition, many of the tissue types included are rare diseases. While many diseases have been profiled in one form or another, this is a comprehensive look across squamous cell carcinomas using proteomics that is similar to the PanCancer analysis of squamous cell from DNA/Methylation/RNA.
Overall this work reports many interesting findings related to the similarities among squamous cell cancers, differences between squamous cell and adenocarcinoma, and differences between common and "rare" squamous cell carcinomas in the context of proteomics. Since this work spans many different squamous tumor types, it is difficult to compare with the many existing findings from the literature. However, Figure 1d provides a nice summary of known targets within the proteomics data.
The paper is primarily a discovery and descriptive paper of many characteristics of squamous cell carcinomas. As such, the work provides both a survey of the proteomics landscape and a resource for other researchers to utilize. The paper provides many complex analyses of the data, providing both the detailed results as supplemental data and figures summarizing the results. Overall, the methods are appropriate although there are a few areas in which there is not enough detail to fully understand the approach taken (see below for specific comments that are likely to be resolved through additional details).

Response:
We appreciate the reviewers for the positive review and valuable comments. We have revised the manuscript according to the comments. The point-to-point responses were as follows.

Q1.
It might be helpful to include p values when describing differences that are qualitatively described within the manuscript. This would allow readers to assess the significance of the observation. For instance (line 100: significant differences across organs) could include a statistical test and associated p value to strengthen the conclusion there are statistically significant differences.
Another example: Stromal vs ESTIMATE "showed a consistent trend" -can this be quantified?  Stromal Ratio (y axis). Spearman correlation. c) Line 267 (revised manuscript): As the statistical analysis was shown in Fig. 3i, we performed Spearman's correlation between FOXO1 and RUNX2 expression. In the revised manuscript, we labeled the R and p value in corresponding results. d) Line 470 and 476 (revised manuscript): For these two sentences "Immunostaining of PRKCE was significantly different among 17 SCCs, and showed an overall high expression in cervical and vagina SCCs" and "In agreement, we noted a high proportion of tumor specific positive SLC27A1 staining in gallbladder and pancreatic SCCs", we showed the immunohistochemistry score for these two markers in pan-SCC cohort, and did Kruskal-Wallis tests (both p < 0.0001) to show the differential expression among 17 SCCs (Figure RL 2).

Q2.
In the case of the overall survival/disease free survival there are cohorts collected at different times having different time to follow up. Therefore, comparing outcomes may be difficult, as censoring may be primarily a late or early event, and may be influenced by date of collection (due to therapy).
There does not appear to be multiple testing correction methods applied to survival p-values. Given the large number of exploratory tests for outcomes, multiple testing correction may be appropriate to consider.
Response: Thank the reviewer very much for giving this critical comment. Due to the complexity of the pan-SCC cohort, we agree that the multiple testing correction is necessary for the survival analysis of the pan-SCC cohort.
The pan-SCC cohort has its characteristics, as it includes 333 patients originating from 17 different organs and ~20 cases per organ. Therefore, we explored the prognostic value of differentially expressed proteins in the pan-SCC cohort. We included age, gender, stage, histology, organ, and protein expression as covariates in the multivariate Cox proportion hazard model and did multiple testing correction (BH adjusted) according to your comments in the revised manuscript. After this strict calculation, we determined molecules with prognostic statistical significance, including Lines 125-126: It would be helpful in showing KRT5 and TP63 to again find a quantification of the assertion that they are "high and ubiquitous". There is some variation in KRT5, at least. It appears that many thyroid cancers have lower KRT5, so it might be helpful to summarize observations with quantification (perhaps by tissue of origin). For instance, how does the variability in AKT1 or KRT5 compare (for instance, what quantile of variability using MAD as done elsewhere)?
Response: Thank the reviewer for these comments. We have calculated the coefficient of variation (CV) and median absolute deviation (MAD) for 333 SCCs separately according to your suggestions ( Table RL1)

Figure RL 3
The protein abundance of SCC diagnostic markers and highly variant genes, the corresponding CV for each marker among 333 SCCs was labeled on the left side.

Q4.
On the basis of the methods and the color gradient in the figure 1d, blue is "missing" data. If this is the case, it may be worth pointing this fact out. FGFR3, PTEN and TP53 may be mutated and thus either not sufficiently translated or translated with alternate forms. Given the (relatively) low expression of those proteins (FGFR3/PTEN) in samples, it may be worth exploring this further.  Table RL 2). Then, we calculated the protein expression level between mutated and wild type in lung SCCs, no statistical difference was found between mutated and wild type SCCs (Wilcoxon rank-sum test, p = 0.55 respectively; Figure RL 4b).
• TP53 mutation and expression in esophageal SCC: The mutation rate of TP53 is 83.3%  From the above analysis, we cannot conclude how mutated FGFR3, PTEN, and TP53 affect the protein expression. One head and neck squamous cell carcinoma study (PMID:33417831) showed that missense mutations in TP53 were associated with increased p53 mRNA and protein abundance, suggesting that specific TP53 mutations might endow oncogenic gain of function to this protein.
Due to the relatively small sample size of a certain type of SCC we tested, a large-scale genomeproteome wide study is needed to better illustrate how gene mutation status affect the protein expression in SCCs.

Figure RL 4 Summarized mutation information associated with protein expression. a TP53, PTEN,
and FGFR3 mutation status in SCC originating from esophagus, lung, and cervix. Comparisons of p53 protein abundance between TP53 mutated and TP53 wild type samples in lung SCC (b) and esophageal SCC (c).
Identification of proteins may relate back to tissue specificity so it would be helpful to know if the total numbers were significantly different across tissues. Specifically, on line 115 "no major differences in the coverage between the 17 SCCs". I would encourage a test of this assertion. Figure 1c suggests there are significant differences. Were samples run by tissue of origin? If they were randomized, then the differences in total protein could be biologically meaningful and worth further investigation. This is particularly true with copy number changes that have been seen in squamous cell cancers.
Response: Thank the reviewer for the comments. It is our fault to state "no major differences in the coverage between 17 SCCs" with no statistics. Our answer is as follows.
a) We performed mass spectrometry profiling randomly, not in the order of organs. c) The thymic SCC and nasopharyngeal SCC are the squamous cell carcinoma with the maximum identification number. As you mentioned, this is probably due to copy number changes in SCCs. We agree that the copy number changes affect protein identification, as tumor samples were identified with more proteins than tumor adjacent normal tissues (PMID: 33417831, 32649877). Copy number changes were frequently detected in SCCs (PMID: 24686850, 31395880). Moreover, the cell density of the thymus and nasopharyngeal SCC is high (as evaluated by HE staining), which may be another reason for the high number of protein identification in these organs.

Q6.
Batch correction: When batch correction was done, there was missing data present. Was the missing data imputed in the batch correction or left missing? That is, the missing data is imputed as a low expression however it is possible that batch correction could increase the level of expression of these missing values? This would presumably introduce an unintended batch effect in that case?
Response: Thank the reviewer for the comments.
• In the input matrix of batch correction, the protein expression value was the original value, the expression value of missing expressed protein was 0, and there was no blank value.
• After batch effect correction, negative values appeared in the matrix. We reset the negative value to 0 (PMID:22257669).
As Figure RL 5a showed, the 293T samples of the SCC and AC cohort were separately clustered before batch effect correction. Batch effect correction was done following the steps described above, and the PCA showed a remarkable similarity between these two batches after batch correction (Figure RL 5b). Therefore, we think this method reduced the introduction of new batch effects. We did the same thing for 7 ACs (4845 proteins).
b) The union set of commonly expressed proteins in SCC (5,130 proteins) and AC (4,845 proteins) contained 5,914 proteins, and the interaction contained 4,061 proteins. A total of 1,069 proteins were only commonly expressed in SCCs, and 46 proteins were only detected in SCCs (not detected in ACs). A total of 784 proteins were only commonly expressed in ACs, and 30 proteins were only detected in ACs (not detected in SCCs).
c) To compare the differences between SCCs and ACs, we removed a total of 76 proteins only expressed in SCCs or ACs, and a total of 5838 proteins were obtained.
In this case, no missing values were in the data matrix, and the batch effect correction was successfully done. Response: We thank the reviewer for the suggestions. Nine TCGA cohorts were originally analyzed by Kaplan-Meier analysis. In the revised version, we did the multiple testing correction (BH adjusted) and modified Fig. 2e, Fig. 2g, supplementary Fig. 6, and supplementary Table   2d with adjusted p values. Fig. 2g shows that kinase and transcription factors played a consistent prognostic role with their downstream targets in SCC or AC, giving an insight that the function of keratinization, glucose metabolism, and extracellular matrix pathway could be consistent or opposite in SCCs or ACs. We want to present this phenomenon, though only in specific tumor types. Now we moved original Fig. 2g to supplementary Fig. 6c. This finding should be explored in future studies.
References were added to the revised paper. Thank the reviewer again for the recommendation. Q9.
Figure 3e: This applies to the PLIN1 finding, but likely applies to a number of other findings. In the case of PLIN1, it appears that the large log2 ratio for rare SCC comes from the fact that most (188/228) common samples did not identify PLIN1. In these cases, it would be helpful to provide any information on presence of peptides associated with PLIN1 in the common cancers. It would be extremely interesting if this protein was only (or substantially) present in rare squamous tumors, potentially due to a deletion in "common" squamous tumors. However due to label-free characteristics it is possible that technical artifacts caused this difference (perhaps fewer proteins overall, run order, peptide interference).
Response: Thank the reviewer for the critical comments. From our result, PLIN1 was highly expressed in rare SCCs and was not detected in the most of common SCCs (188/288) as you mentioned (Figure RL 7a). This is really interesting.
• We firstly checked the database searching result and no PLIN1 peptide were found in common SCC cases with no PLIN1 expression. As you mentioned, we cannot tell the gene status and potentially a deletion in common SCCs.
• Then, we ordered the FISH probe for PLIN1 (Empire Genomics Corp, PLIN1-20-OR) and tested the PLIN1 copy number in ten cases for each SCC. No deletion was found in common SCCs. Interestingly, we detected gene amplification in 3 anal SCCs (3/10, Figure RL 7b).
These results were updated in the revised manuscript.
In this case, we think that the PLIN1 amplification is probably the reason for high expression in rare SCCs. A large-scale study will be needed to explore further the significance of PLIN1 in SCC initiation and progression. Repeated runs with the same samples have high correlation as we expected. Also, this result is consistent with the t-SNE analysis (Fig. 4a), showing that samples from the same organ tended to cluster together. We also included a new Spearman's correlation for all 333 samples in the Supplementary Fig.2c, showing a high correlation within cancer types. Q13.
Supplemental Figure 2c: I do not know what GP Number stands for here, but "cumulative identified proteins" may be more informative.
Response: Thank the reviewer, we agree that "cumulative identified proteins" is more informative than "GP Number" in Supplementary Fig. 2c. In the revised Supplementary Fig. 2e (original Supplementary Fig. 2c), we have changed the "GP Number" to "cumulative identified proteins". Q14.
Supplemental Figure 2g: Is this nomenclature used in the paper?
Response: Thank the reviewer for bringing this to our attention and we only mentioned " Supplementary Fig. 2g" in the submitted manuscript. In the revised version, we have added this nomenclature when referring these protein sets in line 122, 156, 230, and 281.

Q15.
Supplemental Figure 3: This is difficult to understand. For instance, it appears that SPRR1A is high in normal esophagus from the human protein analysis, yet missing (it appears) in the pansquamous cohort. Please clarify this finding.  (Figure RL 9a). It was reported that it had lost expression in esophageal SCC compared to the matched normal tissue samples (Figure RL 9b). As is shown in Figure RL 9c, SPRR1A had lost expression in esophageal SCC. In this case, this proved the tumor purity of the pan-SCC cohort to some degree. We have modified our manuscript to make this clearer. Response: Thank the reviewer for the comments, and it is indeed a little confusing. We labeled Gallbladder SCC, Gallbladder AC, and Breast SCC with red in Supplementary Fig. 5C, as these three cancers grouped in the PCA (Supplementary Fig. 5B). In the revised Supplementary Fig.   5c, we changed the color to black.
Q17. Response: Thank the reviewer for the comments. We have changed the figure legend to "significantly differentially expressed proteins". Q18.
Line 160,161: It is not clear if the term overrepresented is accurate. It appears to be the results of a Wilcoxon test; therefore, it is overexpressed or overabundant.
Response: Yes, we agree and we have changed the "overrepresented" to "overexpressed". Thank the reviewer for the recommendation. Q19.
Line 196-197: It is not clear what data is used to assert that prognosis "were significantly affected by these DEPs".
Response: Thank the reviewer for bringing this to our attention. We have modified this sentence and added "DEPs in ECM, glucose metabolism, and keratinization (Fig. 2e)" to make it clear.

Q20.
Line [194][195]: This appears to be a statement about other work, not specifically results from this cohort although it is not clear.

Response:
Yes, it is a statement of other works and we have moved these findings to Discussion part.
Lines 402-403: it seems likely that the positive rate above 54% is in the HPV+ patient cohort (as is HPV18 6.67), but it is not clear.

Response:
Thank the reviewer for the comments. The total HPV (not type specific) infection rate is ~80%, and HPV16 infection rate is ~54%. This sentence is ambiguous in the submitted manuscript. We have modified this sentence to make it clearer.

Q22.
Lines 410-411: p53 loss occurred in HPV-cases as well. This would be an example where a statistic may provide more information (perhaps the frequency of p53 loss).
Response: Thank the reviewer for the suggestion. We have summarized the p53 loss in HPV positive and negative patients (Table RL 5-7), and labeled the frequency of p53 loss in the revised manuscript. As you mentioned, p53 loss occurred in HPV negative cases as well.
In anogenital SCCs (5 sites,    Line 411: "CDKN2A expression was highly correlated" -please include p value or correlation coefficient. Figure 6b does not include these numbers either. After calculation, we think the statement "CDKN2A expression was highly correlated with HPV infection" is not accurate. We have revised the sentence to "CDKN2A expression was higher in HPV positive perineum SCC than HPV negative perineum SCC". Q24.
Line 418: "We identified 8 patterns of differential pathway ..." -was tissue of origin included in the limma analysis for this? Was this driven by tissue of origin?
Response: Thank the reviewer for the comments. We included the tissue of origin in the limma analysis for eight patterns of the differential pathway, and we think it is more HPV driven patterns.
In this part, we combined all anogenital SCCs, as these organs belonged to the same system and were HPV infection-related. These analyses indicated that HPV16 infection may lead to active inositol phosphate catabolic process and immune evasion, participating in HPV16+ SCC carcinogenesis.
Interestingly, Inositol phosphates were reported promoting HIV-1 assembly and maturation to facilitate viral spread in human CD4+ T cells (PMID: 33476323). Multiple isomers of inositol phosphate were found in Epstein-Barr-virus-transformed (T5-1) B-lymphocytes and may be related with cell transformation or proliferation (PMID: 1660712). We hypothesize Inositol phosphate catabolic process probably related to HPV16 related tumorigenesis.
Q25. refers to the specific tissues selected for staining or in the proteomics cohort?
Response: Thank the reviewer for the comments. Samples in Fig. 7d are all in the proteomics cohort. P16 exhibited positive expression in all HPV+ cervix SCCs and vagina SCCs, and we only selected one representative case to present in Fig. 7d. We also modified the sentence to make it clear.

Q26.
Line 481-482: It would be helpful if there is a test for agreement with the proteomics data and/or the highly expressed in thymus SCC.
Response: Thank the reviewer for the comments. We have added a correlation analysis (Spearman correlation, p =0.015; Figure RL 11) for the immunohistochemistry data and proteomics data.   Additionally, samples were collected from 2001 to present. The manuscript states "all patients provided written consent". This suggests that a general banking protocol was used to collect tissues, distinct from this specific study. It would be helpful to clarify this. In the case of the adenocarcinomas, it is explicitly stated that patients were consented prior to surgery.
Response: Thank the reviewer for bringing this to our attention. A general banking protocol was used to collect tissues and patients were consented prior to surgery. We have moved this statement to the last to make the SCC and AC cohort consistent.

Q30.
Patients with adenocarcinomas were excluded if radiation was used (unclear if this was preoperative). Presumably this was not the case in the squamous cohort, was this evaluated?
Response: Thank the reviewer for bringing this to our attention. The radiation mentioned in the adenocarcinoma part is preoperative. The pan-SCC cohort was also excluded if preoperative radiation was used. We have modified this in the revised manuscript.
Response: Thank the reviewer for the comments. We used the same strict pathological assessment criteria in AC as the SCC. In the revised manuscript, we have added this criterion.

Q32.
Database Searching (640): What version of Uniprot was used? It may be present in a figure, but not stated in the methods. Also, were missed cleavages considered for trypsin identification?
Modifications? More details here would be needed to understand how the searching was done.
Response: Thank the reviewer for the comments and apologize for our negligence. The Uniprot library was downloaded on 2019_07. SwissProt was chosen, including 20,431 human entries. We have modified this in the revised methods. It is true that missed cleavages were considered for trypsin identification, we selected 2 for the missed cleavage. Carbamidomethylation was selected as fixed modification and oxidation was selected as variable modification.
line 759-767: What does "the 20% bottom mad proteins of each organ were combined and duplicates removed" mean? It would seem reasonable to remove the bottom 20% of mad proteins due to lack of variability, but it's not clear what the term "combined" means here (averaged?).
Response: Thank the reviewer for the comments. In this part, we explored the similarities and differences between 17 SCCs.
• Firstly, we'd like to explain why we used "the 20% bottom mad proteins of each organ". To explore the proteomic clustering of 17 SCCs, we screened proteins with consistent and ubiquitous expression in each SCC, meeting the following two criteria. 1, the proteins expressed in at least 1/3 of the samples in specific SCC type; 2, the proteins were sorted according to MAD value. The 20% proteins with the lowest MAD values were selected as the proteins with consistent and ubiquitous expression in each certain SCC. In other words, we think these proteins could represent the molecular features of one certain SCC.
• Next, we "combined" the 17 protein sets, removed the duplicate proteins, and obtained a protein set containing 10,259 proteins. So, "combined" here means that we get a union set of 17 protein sets. In the following steps, we calculated the average expression value of each protein in each SCC, and 1,500 proteins with top MAD values were used to do the hierarchical clustering.
Q34. Response: Thank the reviewer for the comments, and we apologize for the unclear description.
For this part, we'd like to establish a random forest to distinguish these 17 SCCs.
• Firstly, the most consistently expressed proteins were selected within each organ as 1% bottom mad proteins identified in at least 75% cases in each organ. A total of 17 protein lists were then combined, and one protein list containing 1,220 proteins after removing repetitive proteins was obtained (a total of 1,220 proteins were obtained from 333 cases).
• Secondly, all 333 SCC cases were randomly divided into a training set and a validation set, containing 75% and 25% cases, respectively.
• Thirdly, the RandomForest function was used on the training set to calculate the importance (Indicated by "Mean Decrease Accuracy") of each protein (top importance was obtained only from training set). 10-fold cross-validation was used to select the suitable number of top important proteins.
As you mentioned, ordering by mad does not impact the random forest and we have moved this filtering.

Response:
We appreciate your correction and, we have modified this in the revised manuscript.
Response: Thank the reviewer for the suggestion, and this is our negligence. The main type was defined using the minimum cycling threshold in the PCR process when multiple infections happened. We have annotated the main type in the revised manuscript.
line 511: "involved in rare SCC initiation". It would be more accurate to describe them as regulators involved in rare SCC since it is a descriptive finding.
Response: Thank the reviewer for the recommendation. We agree that the 'initiation' is appropriate here, and we have deleted the word.

Reviewer #2 (Remarks to the Author): Expert in SCC proteomics
In the manuscript entitled "Proteomic Landscape of Pan-squamous Cell Carcinomas", Song et al.
provided an interesting panorama of SCC proteomes from 333 patients and 17 sites using mass spectrometry label-free quantitation. By performing several group comparisons (SCC vs AC; common SCC vs rare SCC; HPV+ anogenital SCC vs HPV-anogenital SCC), the authors indicated signatures of differentially abundant proteins that have a prognostic significance and are able to modulate specific pathways. Four groups of SCC types were defined by clustering tumor proteomes and six subtypes of SCCs were determined based on immune enrichment of cell composition. Molecules from these groups may potentially modify biological pathways and are promising as druggable targets. Finally, a machine learning model was built to predict the tumor site of origin based on SCC proteomes and 3 proteins were selected for validation using IHC.
Overall, that is a detailed study that provided insights into the biology of a common cancer type (SCCs) by profiling the tumor proteomes followed by an extensive bioinformatic analysis. The significant number of patients evaluated is certainly a highlight of the study. Since large-scale analysis of tumors available in the literature have focused on DNA/RNA levels, describing the proteome and the potential functional implication is definitely novel and opens up new perspectives for the study of cancer. The description of a protein signature that can possibly define the origin of SCCs is also remarkable and reveled targets with potential to translate to the clinics.
For instance, the authors provided a valuable resource to the scientific communities for further exploration. Considering the methodology, a noteworthy point is that the manuscript presents a deep characterization of SCCs achieved by implementing a strong and robust bioinformatic analysis. The proteomic results indicate a successful workflow of sample preparation, MS run and data analysis (~ 15,000 protein groups quantified for SCCs and ~10,000 for ACs), even though some additional information should be provided to assure the reproducibility and quality of the data.
In summary, the study is of major interest for the oncology and proteomics fields and is in-line with most articles published in Nature Communications. There are therefore some shortcomings that should be addressed.

Response:
We appreciate the reviewers for the positive review and valuable comments. We have revised the manuscript according to the comments. The point-to-point responses were as follows.

Q1.
Line 109: The paper relies on proteomics data from a large group of SCC samples and the quality control of MS runs is certainly one of the main concerns to assure the reliability of bioinformatic  • We plotted the Spearman correlation coefficients for all 333 MS runs as you suggested ( Figure   RL 12a). The median correlation coefficient among these samples was 0.74, and the maximum and minimum values were 0.99 and 0.56, respectively.
• Identified protein number per sample was shown in Figure RL 12b by organ. On average, the SCC proteome had 8,120 protein groups per sample, ranging from a minimum of 6,261 in the thyroid to a maximum of 9,296 in the thymus and 5,648 proteins were present in all 17 SCCs.
• You commented, "if suitable, describing retention times or m/z for trypsin autolysis peaks across samples are suggested approaches." Sorry that the trypsin autolysis peaks were not searched in this work, and we chose the conserved peptides instead. Firstly, we chose to describe the retention times (RTs) of consistently identified peptides for stably expressed proteins, including TP63, EGFR, KRT5, Actin, ELMO2. A total of 65 peptides were identified in all 17 SCCs (not 333 cases), belonging to EGFR, ELMO2, KRT5, and TP63. Secondly, a scatterplot showed the identified peptide frequency and coefficient of variation (CV) of peptide abundance (Figure RL 12c: left). Peptides with a frequency >97% and CV < 0.7 were chosen for display in line chart ( Figure RL 12c: right). As the database search results of Peaks online were only provided RT mean of each SCC, the RT mean of 17 SCCs for four peptides was shown in Figure RL 12c (right).
From these results, we think the works' MS data is of high quality.

Q2.
Normalized protein intensities from pan-SCC and pan-AC were log2 transformed. Even that is not always true, log transformation reduces the skewness of large-scale data and make it more closely to a normal distribution. If that is the case, it would be appropriate to use parametric tests for group comparison instead of the non-parametric analysis described in the manuscript (Wilcoxon rank sum test, Kruskal-Wallis test). Did the authors test the normality assumption of the data? Please comment.
Response: Thank the reviewer for the comments.  3. For instance, SPRR1A, a cross-linked envelope protein of keratinocytes, is an esophageal signature protein (Figure RL 14a). It was reported that it had lost expression in esophageal SCC compared to the matched normal tissue samples (Figure RL 14b). As is shown in Figure RL 14c, SPRR1A had lost expression in esophageal SCC. In this case, this proved the tumor purity of the pan-SCC cohort to some degree. We have modified the result to make this clearer. Response: Thank the reviewer for the comments, and this is our negligence. NES stands for normalized enrichment score, and we have added this explanation in corresponding figure legends.

Q5.
Lines 179 and 325: Besides using their own proteomics data for Cox analysis of DEPs between SCCs and ACs, 9 RNA databases from cancer tissues were used to assume the role of proteins in prognosis. Did the authors evaluate if transcripts are somehow correlated with proteomics data?
Otherwise, the assumptions may not be true.

Response:
Thank the reviewer for this critical comment. Due to the availability of transcriptomics, we validated the proteins using RNA datasets. We did not evaluate the correlation between transcripts and proteomics, as limited corresponding SCC proteomics data was available. Some studies showed positive and significant correlations with the proteomics and mRNA transcripts (PMID: 31585088, 30962452, 31395880). Figure  In addition to the positive correlation between expression of proteins and transcripts, previously published proteomic researches also explored the prognostic value of proteins using survival data of the transcriptome databases, including The Cancer Genome Atlas (TCGA) and Queensland Centre for Medical Genomics (QCMG). One representative study (PMID: 31484774) showed poor prognostic proteins participating in ECM process, such as S100A6 and FN1, which were also captured by our analysis.
In the future studies, we will validate these proteins in SCC proteomics datasets. We talked about the limitation in the revised discussion part.

Figure RL 15
The overall correlation between mRNA and protein data (PMID: 31585088).
Response: Thank the reviewer for bringing this to our attention. In the revised Fig. 5, we have modified this mistake.
Q7. Fig. 5d: Not mentioned in the main text.
Response: Thank the reviewer for the comments. We annotated Fig. 5d in line 326-327 (revised manuscript), but no explanation was provided in the submitted manuscript. In the revised version, we further explain Fig. 5d.
Q8. Fig. 6b: How can the authors explain the high levels of protein pRB in HPV-infected tumors?
Since the silencing of pRb by E7 viral protein produces a rise in p16, the abundance of pRB is not in agreement with what is reported in the literature or in this study. Also, SOX2 is not frequently associated with HPV infection and it would be appropriate to describe the rationale for including this protein in the analysis.
Response: Thank the reviewer for the comments. To answer your question, we firstly compared the RB expression in anogenital vs non-anogenital SCCs. As shown in Figure RL 16a, RB expression is higher in non-anogenital SCCs than anogenital SCCs (Wilcoxon-rank sum test, p < 0.001), which is consistent with previous studies. Moreover, RB showed a significant expression difference among all five anogenital SCCs (Kruskal-Wallis test, p < 0.001) and cervical SCC showed the lowest expression (Figure RL 16b). Spearman correlations between CDKN2A (p16) and RB expressions were then calculated in all anogenital SCCs individually. Interestingly, a negative correlation trend was found in cervical SCC (Spearman correlation; R = -0.255, p = 0.265; Figure RL 16c). So, this cervical SCC data is consistent with previous work.
TP63 and SOX2 were considered as associated with squamous differentiation. Indeed, SOX2 is not frequently associated with HPV infection. We have removed SOX2 from Fig.6b. Thank the reviewer again for this meticulous suggestion. Response: Thank the reviewer for the comments. We added the non-anogenital SCCs in the comparison as you suggested ( Figure RL 17a and 17b).
There was no significant difference between Group1-3 (HPV16 infected cases) and non-anogenital SCC of EZR (Wilcoxon-rank sum test, p = 0.0937), PAWR (Wilcoxon-rank sum test, p = 0.4682), and DUSP3 (Wilcoxon-rank sum test, p = 0.1425) belonging to negative regulation of T cell receptor signaling pathway (Figure RL 17a). As the immune response is ubiquitous in the tissue microenvironment, maybe it is the reason no significant difference between Group 1-3 and nonanogenital SCC.
Interestingly, Inositol phosphates were reported promoting HIV-1 assembly and maturation to facilitate viral spread in human CD4+ T cells (PMID: 33476323). Multiple isomers of inositol phosphate were found in Epstein-Barr-virus-transformed (T5-1) B-lymphocytes and may be related with cell transformation or proliferation (PMID: 1660712). So, we think the Inositol phosphate catabolic process participates in HPV related tumorigenesis. However, due to the small sample size, a large-scale study will be needed to explore this further. Q10.
Line 459: I didn't get the point of why the authors evaluated the 3 markers by IHC. If a signature of 19 proteins was accurately able to discriminate SCC tumor based on their origin (buy the way, this information is not stated in the main text), why the 3 proteins were selected for validation? I understand the validation phase is an important step in large-scale experiments, but analyzing the 3 proteins alone does not make sense in the context of the classifier and did not make the proteomic data stronger.
Response: Thank the reviewer for the comments. We ordered all 19 antibodies to validate the proteomic data, and we successfully bought 16 antibodies. However, 13 antibodies were not getting good staining due to the poor antibody specificity. In this case, we only presented three markers in the study. Now we are ordering clinical grade antibodies and we will validate the other markers in future research when antibodies are available. An explanation had been added in the discussion part.

Q11.
It is not necessary to extensively re-state the key findings in the discussion.
Response: Thank the reviewer for the valuable comments, and we appreciate it. In the revised manuscript, we rewrite the discussion part. On the one hand, reviewers' comments about the discussion part were considered. On the other hand, we discussed the key findings associated with published works instead of just restating these findings.

Q12.
Line 580: References should be provided for WHO classification and TNM system.

Response:
Thank the reviewer for the suggestion. We have found corresponding chapters in the WHO classification and traced the original references. References for SCC classification were added in the revised manuscript (line 597-600, Page 28). TNM systems were referred to the AJCC cancer staging system 8th edition.

Q13.
Fundamental information is missing in the proteomic workflow and it is difficult to judge the reproducibility and quality of the methods employed. Please provide additional information.
Response: Thank the reviewer for the recommendation. We have revised the materials and methods part according to your suggestions. Detailed responses, please see below.
Sample preparation a. Why did the authors add an acetone precipitation step in the FASP protocol? FASP did not take care of contaminant removal? The authors should provide the appropriate references.
Response: Thank the reviewer for the comments. As you mentioned, the FASP did not take care of the contaminant removal very well. In this case, we developed a novel method for FFPE proteomic sample preparation. To decrosslinking and lysing the FFPE samples, we applied a high concentration of detergent (4% SDS) and reducing agent (1mM DTT), which exceeded the normal range in FASP protocol (the number of identified proteins from FFPE samples only using FASP protocol were significantly lower than fresh tissue). Therefore, we added an acetone precipitation step to purify the proteins from the decrosslinking-lysis buffer. Moreover, we used the FASP protocol after the acetone precipitation to dissolve the protein pellets efficiently. The results demonstrated that our protocol was highly repeatable for FFPE proteome profiling with in-depth coverage.
b. How was alkylation and quenching performed?
Response: Thank the reviewer for the comments. We apologize for not describing the alkylation and quenching progress clearly. In our protocol, the homogenized samples were boiled in lysis buffer for decrosslinking and lysis. Acetone precipitation was applied to purify the proteins from the decrosslinking-lysis buffer. Then 8M Urea buffer was used to dissolve the protein pellets. The supernatant was loaded onto an ultrafiltration filter column (10 kD, 500 μl) and centrifuged at 12000 rpm for 15 min. The reduction and alkylation progresses were then carried out, 100 μl of reduction buffer (10 mM DTT, 25 mM NH4HCO3) were loaded on ultrafiltration filter column, incubated for 1 hour at 56℃, and centrifuged at 12000 rpm for 15 min. Then 100 μl of alkylation buffer (55 mM IAA, 25 mM NH4HCO3) were loaded on ultrafiltration filter column, incubated for 45 min in dark at room temperature. The filter was then washed three times by adding 100 μl ammonium bicarbonate (ABC, 50 mM) to the column, followed by centrifugation. The proteins were digested by trypsin. The resulting peptides were loaded on LC-MS.
c. The authors declare that a trypsin-to-protein ratio of 1:50 was used for digestion, but there is no information on protein quantification. How were protein levels determined? The authors also state that "Target on-column load was 200ng total peptide per injection", indicating that peptides were also quantified. Please clarify.
Response: Thanks for the comments. Acetone precipitation was applied to purify the proteins, and 8M Urea buffer was used to dissolve the protein pellets. We measured the protein concentration of the solution using spectrophotometer (NanoDrop, Thermo, USA). The amounts of protein were adjusted to 400 μg. We also measured the peptide concentration of samples before loading to the LC-MS. A total of 200ng peptides was loaded per injection.
d. It is appropriate to present specifications of the trypsin used for digestion. Response: Thank the reviewer for the comments.
• The resolution parameter was set to 50,000 for MS1 and MS2.
• Mass spectra for MS1 and MS2 scans were recorded between 100 and 1700 m/z.
• Ion mobility resolution was set to 0.60-1.60 V·s/cm over a ramp time of 100 ms.
• Data-dependent acquisition (DDA) was performed using 10 PASEF MS/MS scans per cycle with a near 100% duty cycle.
• An active exclusion time of 0.4 min was applied to precursors that reached 20,000 intensity units.
g. Although the identifier is provided in the text, the repository where raw files are deposited is not informed and data could not be accessed.
Response: Sorry for this inconvenience and that is our negligence. In the revised version, we have added the assess link and password.
The accession number for the MS proteomics data reported in this paper is IPX0002831000 Response: Sorry for this inconvenience and thank the reviewer for the comments. the missing values could be significantly reduced. As for the rest missing values after applying for ID transfer, to avoid artificially increasing the false discovery rate, we did not apply other algorithms but a lowest of detection (LOD) strategy for data imputation. We replaced missing values with a certain small number (1/10 of the minimum) to ensure the accuracy of subsequent analysis results. The detailed procedure is as follows. Firstly, we deleted the proteins which were not detected in 2/3 of the samples in each certain SCC. Then, we use the single value (1/10 of the minimum) to replace these missing values. This strategy has been proved to have a robust performance in proteomic data and applied in the previously published studies, such as the proteomic landscape of diffuse-type gastric cancer project (PMID: 29520031) and the early-stage hepatocellular carcinoma project (PMID: 30814741).

Q14.
I was wondering whether the separation of SCCs and ACs in PCA before batch effect correction just reflects their distinct biological characteristics. How can the authors be sure that the separation was a batch effect?
Response: Thank the reviewer for the comments. To confirm the existence of the batch effect, we firstly compared the similarity between the commercial 293T cells of the pan-SCC and pan-AC cohort. As shown in Supplementary Fig. 4a, the 293T samples clustered together within pan-SCC or pan-AC cohort but separated from each other between pan-SCC and pan-AC. From here, we think an actual batch effect existed. Then, we remove the batch effect between pan-SCC and pan-AC using the same method as the 293 samples (Supplementary Fig. 4b).
In this case, we believe the analyses will reflect the differences between SCCs and ACs more accurately.

Q15.
Not clear if IHC was performed in the same cohort as proteomics.
Response: Thank the reviewer for bringing this to our attention. Yes, the IHC was performed using the cases in the proteomics, as we constructed a tissue microarray using the same cohort. We have clarified that the tissue microarrays were constructed using the pan-SCC cohort. Q16.
For HPV grouping, how the authors defined if HPV16 is the main type (group 2) or not (group 3) in multiple infections?
Response: Thank the reviewer for the comments. As the HPV16 infection was tested by RT-PCR in this study, we defined the main type as the one using the minimum cycling threshold. We have revised and added the explanation in line 871.

Reviewer #3 (Remarks to the Author): Expert in SCCs
Understanding the molecular pathways driving histologically similar cancers across anatomic sites may provide new treatment paradigms that historically have been site-specific. Differences in mutational patterns and gene expression profiles between squamous cell carcinomas (SCC) and adenocarcinomas (AC) arising across anatomic sites have been well described using common resources such as TCGA. However, the proteomics landscape of SCC across anatomic sites has not been previously investigated in large numbers of tumors. The current manuscript describes proteomic patterns in 333 treatment-naï ve squamous cell carcinoma (SCC) tissues obtained from 17 organ sites and 69 treatment-naï ve adenocarcinoma (AC) tissues obtained from 7 organ sites from a single university-based hospital in Shanghai, China. Proteomic characterization of tissues was conducted using a mass spectrometry-based approach validated using tissue microarray (TMA) immunohistochemistry (IHC).
Major findings include the elucidation of pathways differentiating SCC from AC (keratinization, glucose metabolism and extracellular matrix), molecules within those pathways associated with disease prognosis, and proteomic clusters/immune subtypes that may represent potential druggable targets. The resulting data repository will serve as a unique and valuable shared resource for investigators to use in the future.
Methods and results are described in great detail, yet there are some key points that should be clarified and/or expanded upon.

Response:
We appreciate the reviewers for the positive review and valuable comments. We have revised the manuscript according to the comments. The point-to-point responses were as follows.

Case selection and classification:
Q1.
The methods state that the cases were randomly selected. How was this accomplished, and what percentage of the total SCC cases treated in the 18-year range do the cases included in the current study represent? Exclusion criteria are presented, yet it's unclear how many patients were excluded for the reasons listed.
Response: Thank the reviewer for pointing out this question, and we apologize for not explaining it clearly. We screened 595 documented SCC patients for 17 organs (Figure RL 18) We also added the patient selection flow chart in Supplementary Fig. 1 (Figure RL 18), and described the detailed information in the methods part as the same time.

Figure RL 18
The quality control and sample filtering standards of sample collection in this cohort.

Q2.
While the overall number of tumors (n=333) is substantial, the numbers of samples available per anatomic site ranged from 10-22 for SCC and 8-12 for AC. Therefore, inferences drawn for specific sites are limited by small sample size. This point should be added to the discussion.
Response: Thank the reviewer for bringing this to our attention and we really appreciate it. This is an important point that should be discussed. The small sample size may limit our findings for specific tumor types, and large-scale studies are needed to validate these findings further. In the revised manuscript, we have added this comment in discussion part. Q3.
The classification of rare versus common tumors is unclear. The authors state that the WHO Classification of Tumors was used, yet there is no reference, and some cancers seem to be misclassified. For example, SCC of the vagina is very rare (i.e. incidence is less than 6 per 100,000), yet it is included here as a common cancer.
Response: Thank the reviewer for the comments, and sorry for this inconvenience.
• The classification of common or rare SCCs in this study was depended on the originated tissue, whether it is squamous epithelium or not.
• Vaginal cancer is a rare gynecologic cancer; however, it was revealed that majority of vaginal cancers reported are SCCs (4 th WHO Classification of tumors of female reproductive organs,

Figure RL 19).
We have made an explanation in our revised manuscript for common and rare SCCs. Also, the references for common or rare SCC classification were inserted in the manuscript.

Figure RL 19
Vaginal SCC epidemiology in 4 th WHO Classification of tumors of female reproductive organs (P211).

Q4.
As the authors point out in the Background, metastatic SCC's (or primaries with elevated metastatic potential) are an important clinical challenge. However, only primary SCC's were included in this case series, and no information was provided on whether or not patients developed metastases during follow-up. This is an important design limitation that should be discussed.
Response: Thank the reviewer for the comments. We pointed out the difficulty in metastatic SCC diagnosis in the background as you mentioned. However, we didn't include metastatic SCCs in this work. The main reason is that we intend to compare the proteome of primary SCCs, and then apply the markers with differentially diagnostic values for metastatic SCCs in the future. We are following up with these patients to acquire the metastases information, and new metastatic SCC cohort is under collecting. In future studies, we'd like to include metastatic SCCs to further validate our findings. In this study, we only compared the characteristics of primary SCCs. This limitation was also discussed in the revised manuscript.

Patient follow-up and survival analysis:
Q5.
No information is provided on the average length of follow-up (and range), as well as whether or not patients were lost to follow-up, and if so, how they were handled in the analysis. It is also not clear that the proportional hazards assumption was assessed.
Response: Thank the reviewer for the comments.
The average length of follow-up is 32 months (3-160 months). For SCCs with a low incidence, 68 patients lost follow-up. These patients were also included in this work but not included in survival analysis.
The multivariate COX proportional hazard model was assessed according to your suggestions.  Please see Supplementary Fig.6 and Fig.4 in the revised manuscript. In future studies, we plan to collect more SCC patients to validate these findings.
In general, it was difficult to follow the results section given the sheer number of figures, figure panels and supplementary materials. The figures were very detailed, as were the supplementary data (often patient-level data files). In many cases, it would have been helpful to create summary tables that allow the reader to directly compare percentages between groups and better ascertain the statistical significance of the observed results. For example, for Figure 1b-it would be helpful to show the information in tabular form so that percentages of samples across anatomic sites could be more directly compared with respect to tissue characteristics such as stromal score and keratinization; statistically significant testing could be used to determine which differences are most likely to be real and not due to chance. The raw data are included in supplementary  Fig.6c to show HPV16 related five Groups distribution in 5 anogenital SCCs. Q8.
Regarding the HPV results, it is difficult to glean from Figure 6e whether the protein patterns depicted are specific to HPV 16. It would be helpful to present HPV type-specific prevalence by tumor type in tabular form, and then present the percentages of each of the five HPV groups defined in Fig 6c that express proteins corresponding to the different pathway groups of interest defined in 6d. Were HPV16 E6 and E7 proteins detected in any of the SCC samples?
Response: Thank the reviewer for giving this recommendation. We added two tables according to your suggestions. Table RL 9 showed HPV-type specific prevalence in 5 anogenital SCCs, and Table RL 10 referred to five groups distribution in 5 anogenital SCCs.
• Then we further explored the protein expression in Inositol phosphate catabolic process as you suggested (Figure RL 20c). INPP1 showed a high expression in Group1 compared to other groups (Figure RL 20c). Interestingly, Inositol phosphates were reported promoting HIV-1 assembly and maturation to facilitate viral spread in human CD4+ T cells (PMID: 33476323).
Multiple isomers of inositol phosphate were found in Epstein-Barr-virus-transformed (T5-1) B-lymphocytes and may be related with cell transformation or proliferation (PMID: 1660712).
We hypothesize Inositol phosphate catabolic process probably related to HPV related tumorigenesis. However, due to the small sample size, a large-scale study will be needed to explore this further.
To answer the question "whether HPV16 E6 and E7 proteins were detected in any of the SCC samples", we did two things.
• On the one hand, we ordered the antibodies for HPV16 E6 (Abcam, Clone: ab70) and E7 (ProMab, Clone: 6F3) and immunohistochemistry (IHC) was performed on all pan-SCC samples. After evaluating IHC staining, we believe that the antibodies with high background and poor specificity, as HPV negative samples were with positive staining (dilution 1:10000 for both antibodies).
Response: Thank the reviewer for this valuable comment. We think that fewer proteins were predictive of survival in the pan-SCC cohort, probably due to the cohort's complexity and the survival analysis method we used is really strict.
• The pan-SCC cohort included a total of 333 patients from 17 different SCCs, while each of the TCGA dataset is just one tumor type.
• The survival analysis method for pan-SCC we used in this study is the multivariate Cox proportional hazards model. We included age, gender, stage, histology, organ, and protein expression as covariates in the revised analysis as you suggested in Q6 and did multiple testing correction (BH adjusted) per reviewer 1's comments. In this case, the numbers of proteins with prognostic values becomes less. The survival analysis method for TCGA cohorts is Kaplan-Meier curve with log-rank test.

REVIEWERS' COMMENTS
Reviewer #1 (Remarks to the Author): This paper revision addresses the concerns and questions. We appreciate the significant effort and additional analysis that was performed in support of this revision.
Several issues still should be addressed: -Batch effect analysis: Line 741, "interaction" should be "intersection".
-Please clarify in the text line 97 "68 patients lost follow-up". Does this mean that these patients did not have any information? The supplemental table suggests these patients are missing information. This is different from lost to follow-up which assumes information about the patient until a specific time point (at which the patient outcome is censored). As there are censored patients in the KM plots, I believe these patients have no outcome information. If so, "68 patients with no outcome information" or similar would be more appropriate.
-Line 462-464: Can you verify that the accuracy was from the training set (the manuscript states this)?
Although it is stated that the data was split into 75/25 there is no accuracy reported for the validation set? Figure 7b shows 100% sensitivity/specificity (as reported from the training set) but the counts of tumors (Tumors, No) suggests that results may be based on validation set. Accuracy on the training set is not as informative as it is expected to be optimistically biased, but it is stated this way in the manuscript. However, if the intent is to report on the validation set accuracy this was not described.